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ABSTRACT 

We assemble the equations of general relativistic magnetohydrodynamics (MHD) in 3 -I- 1 form. These 
consist of the complete coupled set of Maxwell equations for the electromagnetic field, Einstein’s equa¬ 
tions for the gravitational field, and the equations of relativistic MHD for a perfectly conducting ideal 
gas. The adopted form of the equations is suitable for evolving numerically a relativistic MHD fluid in 
a dynamical spacetime characterized by a strong gravitational field. 


1. INTRODUCTION 

Magnetic fields play a crucial role in determining the 
evolution of many relativistic objects. In any highly con¬ 
ducting astrophysical plasma, a frozen-in magnetic field 
can be amplified appreciably by gas compression or shear. 
Even when an initial seed field is weak, the field can 
grow in the course of time to significantly influence the 
gas dynamical behavior of the system. If, in addition, 
the gravitational held is strong and dynamical, magnetic 
helds can even affect the entire geometry of spacetime, 
according to general relativity. In this situation, terms 
involving magnetic and electric helds are important not 
only as electromagnetic forces acting on the matter in the 
equations of relativistic hydrodynamics but also as stress- 
energy sources governing the metric in Einstein’s gravita¬ 
tional held equations. 

In this paper (hereafter Paper I) we assemble the 
complete set of MaxwelUEinstein-magnetohydrodynamic 
equations that determine the self-consistent evolution of 
a relativistic, ideal MHD huid in a dynamical spacetime. 
Our goal is to set down a formulation of the equations that 
is suitable for numerical integration in full 3-1-1 dimen¬ 
sions. Subsets of these equations have appeared elsewhere, 
but we recompile the complete set here for convenience 
and future reference. We reconcile several seemingly dif¬ 
ferent, but equivalent, forms of the equations that have 
appeared in the existing literature. We also correct some 
errors in previously published results. In a companion pa¬ 
per (Baumgarte & Shapiro, 2002, hereafter Paper H) we 
use these general relativistic MHD equations to follow the 
gravitational collapse of a magnetized star to a black hole. 

We are motivated in part by the growing list of im¬ 
portant, unsolved problems which involve hydromagnetic 
effects in strong-held dynamical spacetimes. The hnal fate 
of many of these astrophysical systems, and their distin¬ 
guishing observational signatures, hinge on the role that 
magnetic helds may play during the evolution. Some 
of these systems are promising sources of gravitational 
radiation for detection by laser interferometers now un¬ 
der design and construction, like LIGO, VIRGO, TAMA, 
GEO and LISA. Others may be responsible for gamma- 
ray bursts. Recent examples of astrophysical scenarios in¬ 
volving strong-held dynamical spacetimes in which MHD 


effects may play a decisive role include the following: 

• The merger of binary neutron stars. The merger can 
lead to the formation of a hypermassive star supported by 
differential rotation (Baumgarte, Shapiro & Shibata 2000; 
Shibata & Uryu 2000). While such a star may be dynam¬ 
ically stable against gravitational collapse and bar forma¬ 
tion, the radial stabilization due to differential rotation is 
likely to be temporary. Magnetic braking and viscosity, 
driven by differential rotation, combine to drive the star 
to uniform rotation, even if the seed magnetic held and the 
viscosity are small (Shapiro 2000). This process inevitably 
leads to delayed collapse, which will be accompanied by a 
delayed gravitational wave burst. 

• Core collapse in a supernova. Core collapse may again 
induce diherential rotation, even if the rotation of the pro¬ 
genitor at the onset of collapse is only moderately rapid 
and almost uniform (see, e.g., Zwerger & Muller 1997; 
Rampp, Muller & Ruffert 1998, and references therein). 
Differential rotation can wind up a frozen-in magnetic held 
to high values, at which point it may provide a signihcant 
source of stress. Hypermassive neutron stars may form and 
survive until the helds are weakened by magnetic braking 
or other instabilities. 

• The generation of gamma-ray bursts (GRBs). Typical 
models for CRB formation involve the collapse of rotating 
massive stars to a black hole (MacFadyen & Woosley 1999) 
or the merger of binary neutron stars (Narayan, Paczyn- 
ski & Piran 1992) or the tidal disruption of a neutron 
star by a black hole (Ruffert & Janka 1999). In current 
scenarios, the burst is powered by the extraction of rota¬ 
tional energy from the neutron star or black hole, or from 
the remnant disk material formed about the black hole 
(Vlahkis & Konigl 2001). Strong magnetic helds provide 
the likely mechanism for extracting this energy on the re¬ 
quired timescale and driving collimated CRB outhows in 
the form of relativistic jets (Meszaros & Rees 1997; Sari, 
Piran & Halpern 1999; Piran 2002). Even if the initial 
magnetic helds are weak, they can be amplihed to the re¬ 
quired values by diherential motions or dynamo action. 

• Supermassive star (SMS) collapse. SMSs may form 
in the early universe and their catastrophic collapse may 
provide the origin of supermassive black holes (SMBHs) 
observed in galaxies and quasars (see Rees 1984 and Baum- 
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garte & Shapiro 1999a for discussion and references) If a 
SMS is uniformly rotating, cooling and secular contraction 
will ultimately trigger its coherent dynamical collapse to a 
SMBH, giving rise to a burst of gravitational waves (Saijo 
et al. 2002; Shibata & Shapiro 2002). If a SMS is differ¬ 
entially rotating, cooling and contraction will instead lead 
to the unstable formation of bars or spiral arms prior to 
collapse, and the production of quasi-periodic waves (New 
& Shapiro 2001a,b). Magnetic fields and turbulence pro¬ 
vide the principle mechanisms that can damp differential 
rotation in such stars (Zel’dovich & Novikov 1971; Shapiro 
2000) and thus determine their ultimate fate. 

• The r-mode instability in rotating neutron stars. 
This instability has recently been proposed as a possi¬ 
ble mechanism for limiting the angular velocities in neu¬ 
tron stars and producing observable quasi-periodic grav¬ 
itational waves (Andersson 1998; Friedman & Morsink 
1998; Andersson, Kokkotas & Stergioulas 1999; Lindblom, 
Owen & Morsink 1998). However, preliminary calculations 
(Rezzolla et al. 2000, 2001a,b and references therein) sug¬ 
gest that if the stellar magnetic field is strong enough, 
r-mode oscillations will not occur. Even if the initial field 
is weak, fluid motions produced by these oscillations may 
amplify the magnetic field and eventually distort or sup¬ 
press the r-modes altogether. [R-modes may also by sup¬ 
pressed by non-linear mode coupling (Arras et al. 2002; 
Schenk et al. 2002).] 

This paper is partitioned as follows: In Sections 2 and 3 
we review Einstein’s field equations and Maxwell’s equa¬ 
tions in 3 -|-1 form. In Section 4 we discuss the approxima¬ 
tion of ideal magnetohydrodynamics and in Section 5 the 
equations of general relativistic hydrodynamics. In Sec¬ 
tion 6 we then develop the equations of general relativistic 
MHD. We derive the MHD source terms that appear in 
Einstein equations in Section 7. We compare our results 
with those of previous treatments in Section 8. Finally, we 
briefly summarize our analysis in Section 9. 

We adopt geometrized units throughout, setting G = 
1 = c, where G is the gravitation constant and c is the 
speed of light. 

2. EINSTEIN’S FIELD EQUATIONS IN 3-|-l FORM 

The spacetime geometry (i.e. metric) is determined by 
integrating Einstein’s field equations of general relativity. 
Most algorithms for performing this integration numer¬ 
ically are based on a 3 -I- 1 decomposition of Einstein’s 
equations, which is ideally suited for solving the general 
initial value problem. Below we briefly summarize the key 
equations that result from recasting Einstein’s equations 
in 3 -|- 1 form. More detailed discussions may be found, 
for example, in Misner, Thorne & Wheeler (1973), York 
(1979), Evans (1984) and references therein. 

In a 3 -I- 1 decomposition of Einstein’s field equations, 
the four-dimensional spacetime M is foliated into a family 
of non-intersecting spacelike three-surfaces S, which arise, 
at least locally, as level surfaces of a scalar time function t. 
The spatial metric jab on the three-dimensional hypersur¬ 
faces E is induced by the spacetime metric gab according 
to 

Jab = gab + naUb, (2-1) 

where n“ is the unit normal vector Ua = aVat to the slices. 
Here the normalization factor a is called the lapse func¬ 


tion. The time vector is constructed so that it is dual 
to the foliation 1-forms Vat: 

= ( 2 - 2 ) 

where the shift vector /3“ is spatial, i.e., na/3“ = 0, but 
otherwise arbitrary. In a coordinate system that is aligned 
with t“ and E, the components of n“ are 

= (—a, 0, 0, 0) and n“ = a“^(l, —/?*). (2-3) 

We adopt the convention that roman indices a, b,c,d,... 
denote spacetime components, while denote 

spatial components. 

The spacetime metric can now be written in the ADM 
form (Arnowitt, Deser & Misner 1962) 

ds^ = —a^dt^ + Jij{dx^ + P^dt){dx^ + P^dt). (2-4) 


The lapse function a determines by how much proper time 
advances along the normal vector from one time slice to 
the next, and the shift vector /3* determines by how much 
spatial coordinates are shifted on the new slice. The lapse 
function and three components of the shift vector consti¬ 
tute gauge potentials that may be freely specified. To¬ 
gether, a and /3* thus embody the four degrees of coordi¬ 
nate freedom inherent of general relativity. 

Einstein’s equations 

Gab = STrTab, (2-5) 

where Gab is the Einstein tensor associated with gab and 
Tab is the stress-energy tensor, can be projected both along 
the normal direction n“ and into the spatial slice E. The 
spatial projection yields two constraint equations, which 
constrain the fields within each slice E; they contain at 
most one time derivative of the spatial metric. The pro¬ 
jection along the normal vector yields an evolution equa¬ 
tion that describes how the fields propagate from one slice 
to the next; it contains second-order time derivatives of 
the spatial metric. The constraint equations consist of the 
Hamiltonian constraint 

R + K^ - = Wnp (2-6) 

and the momentum constraint 

- j^^K) = SttS'*, (2-7) 

and the evolution equation is 

dtK,, = -D,D,a + a{R^j-2K,kKk+KK,,) 

-8TTa{Sij - ^jij{S - p)) + CpKij. 

... ( 2 - 8 ) 

Here Kij is the extrinsic curvature, and its definition in 
terms of the time derivative of the spatial metric is usu¬ 
ally considered the second evolution equation 

dtjij = -2aKij + Cpjij. (2-9) 

Here D^, Rij and R = j^^ Rij are the covariant deriva¬ 
tive, Ricci tensor and scalar curvature associated with jij, 
while K = j'^^Kij is the trace of the extrinsic curvature. 
The symbol C denotes a Lie derivative. The matter and 
nongravitational field sources p, Si and S^J are the pro¬ 
jections of the stress-energy tensor into n“ and E and are 
given by 

p = nanbT°’^ 

5 . = -J^aUbT^^ 

Sij = JiaJjbT^’^- 


( 2 - 10 ) 

( 2 - 11 ) 

( 2 - 12 ) 
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The quantity p is the total mass-energy density as mea¬ 
sured by a normal observer, Si is the momentum density 
and Sij is the stress. Finally, S is defined as the trace of 

S = (2-13) 

We remark that if the constraint equations are satished on 
an initial time slice E, the evolution equations guarantee 
that the constraints will be satisfied on all subsequent time 
slices. 

Equations (2-6) through (2-9) are commonly referred 
to as the ADM equations (Arnowitt, Deser & Misner, 
1962). Numerical implementations of these equations have 
revealed that their numerical stability can be improved 
dramatically by bringing them into a slightly different 
form. One such modihcation, now commonly referred to 
as “BSSN”, is based on Shibata & Nakamura (1995) and 
Baumgarte & Shapiro (1999b). This system is widely used, 
and its enhanced stability properties have been analyzed 
by several authors (e.g. Alcubierre et ah, 2000; see Knapp, 
Walker & Baumgarte 2002 for an electromagnetic anal¬ 
ogy). Alternatively, several authors have experimented 
with hyperbolic formulations of Einstein’s equations (e.g. 
Anderson & York 1999). We refer the reader to these pa¬ 
pers for further details and references. 

3. maxwell’s equations 

We decompose the Faraday tensor as 

F<^b ^ 

so that and 5“ are the electric and magnetic fields 
observed by a normal observer n“. Both Helds are purely 
spatial, whereby 

= 0 and = 0, (3-2) 

and the three-dimensional Levi-Civita symbol Cabc is de¬ 
fined by 

^abc ^ ^ u'^edabc- ( 3 - 3 ) 

Note that e®*"’ is non-zero only for spatial indices, while 
Cabc may be non-vanishing even if one index is timelike 
(see equation (3-19) below). We also decompose the elec¬ 
tromagnetic current four-vector according to 

= nVe + (3-4) 

where pe and J“ are the charge density and current as ob¬ 

served by a normal observer n“. Note that </“ is purely 
spatial, J“na = 0 . 

With these definitions. Maxwell’s equations, 

= dTT (3-5) 

and 

V[aF6c]=0, (3-6) 

where V is the four-dimensional covariant derivative op¬ 
erator associated with gab^ can be brought into the 3-1-1 
form 

AE* = dTTpe (3-7) 

dtE^ = F^^DjiaBk) - 47raJ* -fi aKE^ + CpE^ (3-8) 

DiB^ = 0 (3-9) 

dtB^ = -F^^DjiaEk) + aKB^ + CpB^ (3-10) 

(see, e.g., Thorne & MacDonald 1982). The charge con¬ 
servation equation, 

(3-11) 


which is implied by (3-5), becomes 

dtPe = -Di{aF) + aKpe + Cppe- (3-12) 

The special relativistic Maxwell’s equations can be recov¬ 
ered very easily by evaluating the above equations for a 
Minkowski spacetime with jij = fij, where fij is the fiat 
spatial metric in an arbitrary coordinate system, a = 1 , 
a: = 0 and /3* = 0 . 

It is convenient to introduce a four-vector potential Aa, 
which can be decomposed into 

Aa = ^ria + Aa, (3-13) 

where Aa is purely spatial, Aan“ = 0. Inserting (2-3), this 

implies 

At = /3M, (3-14) 

while A* = 0, as for any spatial vector. In terms of the 
vector potential, the Faraday tensor can be written 

Bab — Ab^a Aa^b — ^aBb Tl^Ea F (^abcB (3-15) 
Contracting this equation with yields 

e“'’=(A.a - Aa,b) = e'^^’FabdB'^ = 2B- (3-16) 

or 

B^ = B^^Ak,j. (3-17) 

Note that with this identification, the magnetic field R* 
automatically satisfies the constraint equation (3-9). 

It is possible, and often convenient (see Paper II), to 
re-write Maxwell’s equations completely in terms of Ei 
and Ai, thereby eliminating A. Evaluating (3-15) for 


the components a = t and b = i with At = Ai and 
At = —q; 4> -I- yields 

dtA, = -aE, + eujB^ - (a$ - AA,)„ . (3-18) 

Using the definition (3-3), we can rewrite 

^tij — H ^dtij — bX (3 Cktij — bX /? Ctikj 

= -P^n’^edtkj = -f3’"e,kj, (3-19) 

so that 

dtA, = -aE, + CijkP^B'^ - (a$ - AA,)„ . (3-20) 

With (3-17), eijkP^B’^ can be expressed in terms of Ai as 

eijkfi^B'^ = 

= (^',< 5 ™ - 5^A^^i 

= l3^Aj^,-P^Ai^j. (3-21) 

Inserting this into (3-20) yields 

dtAi = -aEi - 9,(a$) + CpAi. (3-22) 


In equation (3-8), the magnetic field A can be eliminated 
similarly: 

B^'^D,{aBk) = B^^D,{aekimD^A^) 

= {S\S^m - S\S\)D,{aD^A^) 

= Dj{aD*A^)-Dj{aD^A^). (3-23) 
Inserting this into (3-8) yields 

a,A = Dj{aD^A^)-Dj{aD^A^)-ATTar + aKE* + CpE^ 

(3-24) 

Equations (3-22) and (3-24) form a system of equations 
for A* and Ai alone. In the special relativistic limit, they 
again reduce to familiar expressions. We also note that 
in terms of partial derivatives, equation (3-24) can be ex¬ 
panded to yield 

a*A = 7 - 1/2 

-47raA -k aAA -k CpE\ (3-25) 

where 7 is the determinant of the spatial metric jij. This 
form of the electric field evolution equation will be useful 
for applications in Paper 11. 


Va = 0 
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4. IDEAL MAGNETOHYDRODYNAMICS APPROXIMATION 


5. GENERAL RELATIVISTIC HYDRODYNAMICS 


Ohm’s law can be written (see, e.g., Jackson, 1999) 

Ja - PeUa = CjFabU^, (4-1) 

where cr is the electrical conductivity and pe = —J°'Ua 
is the charge density as seen by an observer co-moving 
with the fluid four-velocity (in contrast to pe, which 

was defined as the charge density as observed by a normal 
observer n“). 

A 3 -I- 1 decomposition of Ohm’s law can be derived by 
contracting (4-1) with n“ and The former yields 

Wpe =Pe- aUaE\ (4-2) 

where we have defined W as the Lorentz factor between 
normal and fluid observers 

W = —naU°' = au*. (4-3) 

Projecting (4-1) into E, or, equivalently, evaluating the 
spatial components a = z of (4-1), yields 

Ji-pe.Ui = aFiaU°‘= a{FitU* + FijU^) 

= a{aE^u* + eukB’^u* + 

= a (WE, + B’^u* + B’^u^)) 

= a (WE, + eiju(v^ + l3i)B^u^)) . 

Here we have defined 


(4-4) 


(4-5) 

F 

and have used (3-19) to relate tnk to e^fc, giving rise to 
the shift term in (4-4) (the shift term is missing in some 
previous treatments, see Section 8). 

Dividing Ohm’s law (4-1) by a and allowing cr —> oo 
yields the perfect conductivity condition 

Eabv!’ = 0. (4-6) 

According to (4-2) and (4-4), this result is equivalent to 
the condition that the electric field vanish in the fluid rest 
frame 

= 0, (4-7) 


or 

ocEi = -e,jk(F +13^)B\ (4-8) 

which is often called the ideal MHD relation. When eval¬ 
uated in a Minkowski spacetime, the last equation reduces 
to the familiar expression Et = —eij^v^B^ or E = — v x B. 

We can now evaluate Faraday’s law (3-10) under the as¬ 
sumption of perfect conductivity. Taking the trace of (2-9) 
yields 

aK =+ Di(i\ (4-9) 

The above expression can be combined with the Lie deriva¬ 
tive CfjB^ to give 

aKB^ + Cf)B^ = Dj{(3^B^ - f3^B^) - B^dt (4-10) 

where we have used (3-9). Inserting (4-10) together with 
the ideal MHD equation (4-8) into Faraday’s law (3-10) 
reveals that all the shift terms cancel, leaving 

^dt(i^/^B^) = D,(FBi - FB^). (4-11) 

It is convenient to introduce the magnetic vector density 

(4-12) 

in terms of which equation (4-11) and (3-9) reduce to the 
particularly simple forms 

dtB^ = dj(v^B^ - v^B^), (4-13) 

and 

d,B^ = 0. (4-14) 

In Section 8 we compare with previous treatments and 
correct errors in some previously published equations. 


For a perfect fluid, the stress-energy tensor Tg can be 
written 

= + (5-1) 

Here po is the rest-mass density as observed by an observer 
co-moving with the fluid zz®, P is the pressure, and h is 
the specific enthalpy 

h = l-|-e-|-P/ pq, (5-2) 

where e is the specific internal energy density. 

In the absence of electromagnetic fields, the equations 
of motion for the fluid can be derived from the local con¬ 
servation of energy-momentum, 

VftTfltg = 0, (5-3) 

and the conservation of baryons, 

Va(pou“) = 0. (5-4) 

The resulting equations can be cast in various forms, de¬ 
pending on how the the primitive fluid variables are chosen 
(see, e.g., Font 2000 for a recent review). The most fre¬ 
quently adopted relativistic formalism was originally de¬ 
veloped by Wilson (1972; see also Hawley, Smarr & Wil¬ 
son, 1984), who defined a rest-mass density variable 

D = poW, (5-5) 

an internal energy density variable 

E = poeW, (5-6) 

and a momentum variable 

Sa = PohWua = (D + E + PW)Ua. (5-7) 

Note that the spatial vector Si defined above is the fluid 
contribution to the source term Si appearing in Einstein’s 
field equations (cf., equation (2-11)). In terms of these 
variables, the equation of continuity becomes 

^((yi/^P) -h dj(j^^‘^Dv^) = 0. (5-8) 

Contracting (5-3) with F yields the energy equation 

dt(j^/^E)+d,(j^/^Ev^) = 

- P (dt(-i^/^W) + di(-i^/^WF)) , (5-9) 

while the spatial components of (5-3) yield the Euler equa¬ 
tion 

dt(i^/^S,) + dj(^^/^S,F) = 

(5-10) 

For gamma-law equations of state 

P=(r-l)poe, (5-11) 

the right hand side of the energy equation (5-9) can be 
eliminated to yield 

dt(l^/^E,) + d,(-i^/^E,F) = 0, (5-12) 
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where we have introduced the energy variable E* dehned 
as 


E* = {poef/^W (5-13) 


may exchange energy and momentum. In particular, one 
hnds 


VfeTfltd = (6-2) 


(see Shibata 1999, who also absorbed the determinant 
7 ^/^ into the definition of the fluid variables; see also Shi¬ 
bata, Baumgarte & Shapiro 1998). This simplification has 
great computational advantages, since the time derivatives 
on the right hand side of (5-9) are difficult to handle in 
strongly relativistic fluid flow (compare Norman & Win¬ 
kler 1986). 

For given values of Uj, W can be found from the nor¬ 
malization relation UaU°' = — 1 , 

W = au* = (l -I- , (5-14) 

and w* from 

= af^UjlW - /3\ (5-15) 

Finite difference implementations of the above equations 
must be adapted to handle the appearance of shock dis¬ 
continuities. Two strategies are commonly adopted. 

The more traditional approach is to add an artificial 
viscosity term to the equations (von Neuman & Richt- 
myer 1950). Typically, the artificial viscosity term Q is 
non-zero only where the fluid is compressed and is added 
to the pressure on the right hand sides of both the energy 
equation (5-9) and the Euler equation (5-10). The arti¬ 
ficial viscosity spreads the shock discontinuity over sev¬ 
eral grid zones. For shocks occurring in Newtonian fluids 
with modest Mach numbers, artificial viscosity generates 
the Rankine-Hugoniot jump conditions to reasonable accu¬ 
racy. Artificial viscosity has also been used successfully in 
relativistic applications (e.g. Wilson 1972; Shibata 1999), 
but it leads to less satisfactory results for highly relativistic 
flows or high Mach numbers (Norman & Winkler 1986). 

An alternative approach to handling shocks is a high res¬ 
olution shock capturing (HRSC) scheme (see, e.g., Marti 
& Miiller 1999 for a recent review). In such a scheme, one 
treats all fluid variables as constant in each grid cell. The 
discontinuous fluid variables at the grid interfaces serve as 
initial conditions for a local Riemann shock tube problem, 
which can be solved either exactly or approximately. Al¬ 
lowing for discontinuities, including shocks, lies at the core 
of these schemes, and does not require any additional ar¬ 
tificial viscosity. Constructing Riemann solvers for HRSC 
requires knowledge of the local characteristic structure of 
the equations to be solved. This has motived the develop¬ 
ment of several flux-conservative hydrodynamics schemes, 
which do not contain any derivatives of the fluid vari¬ 
ables in the source terms, and for which this characteristic 
structure can be determined (see, e.g., Font 2000; Font et 
al. 2002 ). 

6. GENERAL RELATIVISTIC MAGNETOHYDRODYNAMICS 

To derive the equations of general relativistic magne¬ 
tohydrodynamics, we now add the electromagnetic stress- 
energy tensor to the fluid stress-energy tensor 

Tab = TCd + T:^- ( 6 - 1 ) 

Local conservation of energy-momentum demands that the 
divergence of the sum Tab vanish. The divergence of the 
individual stress-energy tensors T^^id ^^4 do not van¬ 
ish in general, since the fluid and electromagnetic fields 


(see equation (5.40) in Misner, Thorne & Wheeler 1973). 
The right hand side of equation (6-2) now includes the 
Lorentz force in the equations of relativistic hydrodynam¬ 
ics. Note that the baryon conservation equations (5-4) and 
(5-8) remain unchanged. 

In the energy equation (5-9), which was derived from 
the contraction UhVaT®**, the addition of the Lorentz force 
yields 

= (6-3) 

- P 

For a gamma-law equation of state the above equation may 
be written in terms of F* according to 


at(7^/^F*) -I- = 


- UaF’^'^Jb 





(6-4) 


which now takes the place of (5-12) in general. The Euler 
equation (5-10) now becomes 

5t(7'/'5,) + a,(7'/"5,T^')= (6-5) 

- (^iP + -k OL-i^^'^FiaJ°‘. 

In the case of ideal MHD the new terms on the right 
hand sides of the energy equations (6-3) and (6-4) vanish 
because of (4-6). This result is understandable, since it 
corresponds to the absence of Joule heating in the limit of 
infinite conductivity. 

We now proceed to determine the Lorentz force 
in the Euler equation (6-5). Since 77“ is not known a priori, 
we will first use (3-5) to express Ja in terms of the elec¬ 
tromagnetic fields. Note that Ja is the current four-vector 
as opposed to its spatial projection ,T. We could express 
E immediately using the spatial Maxwell equation (3-8). 
Instead we need to derive the four-dimensional equivalent 
of (3-8) to express Ja, which we then can contract with 
the Faraday tensor (3-1) to obtain the Lorentz force. 

The divergence of the Faraday tensor in (6-2) is 

VhF“*' = n“V6F*' -h F'’Vhn“ - n^VbE’^ - E’^V^n^ 

+Vhe“*'“R,. (6-6) 

We will now decompose the four-dimensional derivatives 
Va in each term above. 

The Lie derivative of along arE is 

= an*'VhF“ - aE'^VbrP - E’’n‘^Vba, (6-7) 
or, with (2-2), 

- (dt - Cb)E'^ = n‘’\/bE^ - E^\/bn^ - E'^n^Vb In a. (6-8) 
a 

The four-dimensional divergence VaE°' can be expressed 
in terms of the three-dimensional divergence DiE^ 

DaE<^ = -la^VbE- = {gj> + nan^)VbE- 

= VaF“-F“D,lna, (6-9) 

where we have used naF“ = 0 and 

n^VbUa = Oa = Da In a. 


( 6 - 10 ) 
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Here Qa is the four-acceleration of a normal observer. Since 
the extrinsic curvature Kab can be written 

Kab = -VaUb - riaab, (6-11) 

the divergence of n“ satisfies 

Van“ = -K. (6-12) 

Inserting these expressions into (6-6), we now find the in¬ 
termediate result 

= n^DbE^’ + KE^ - -(5* - Cb)E'^ 
a 

(6-13) 

where we have also used E°‘Da In a = E°‘Va In «■ 

The term involving Ba in (6-13) can be rewritten as 

bB, + B.VtUd) 

= e°^^‘^‘’-{ndDbBc - BcUbad) 

= a~^e°'^'^'^{andDbBc + ndBcDta) 

= a-h^'^^DbiaBc), (6-14) 

where we have used equations (3-3), (6-10) and (6-11). 
Inserting this expression into (6-13) now yields 

47raJ"“ = aVhF“*' = -{dt - Cb)E^ + e‘^^'^Db{aB^) 

an’^DbE^ + aKE‘^. (6-15) 

Not surprisingly, this is the four-dimensional version of 
equation (3-8), which can be found by taking the spatial 
projection of (6-15). 

The next step is to contract (6-15) with the Faraday ten¬ 
sor (3-1). Using noA“ = 0, noe“^° = 0 and UaCanE^ = 0 
several terms cancel, and one finds 

^^aJ^Fab = aEaDbE^ (6-16) 

+ Ua {-Ebidt -Cb- aK)E^ + Ebe^^<^D,{aBd)) 

- B-ecab {{dt -Cb- aK)E'^ - e^^^Dd{aBe)) . 

This expression can now be inserted into the Euler equa¬ 
tion (6-5). For spatial components rii = 0, so that the sec¬ 
ond line in (6-16) vanishes. The source term a^^l'^FiaJ°' 
can then be rewritten 

1 /2 

= ^( - B^t,,k{dt -Cb- aK)E^ 

+ aE,{DjE^) + BWj{aBi) - BW,{aBj)^ . (6-17) 

If desired, the covariant derivatives in the last two terms 
can be converted into partial derivatives, which finally 
yields the Euler equation 

- {d.P + 

- {eMdtE'^ - f3'diE^ + E^difi’^ - aKE>^) 

+ di{aBj) - dj{aBi)^ , (6-18) 

where we have expanded the Lie derivative of AL Note 
that in this equation the electric field terms enter with the 
opposite sign from those in the corresponding equation 
(3.2) of Sloan & Smarr (1985, hereafter SS), who further 
assume (3 = 0 = K. 


For numerical implementations, the most challenging 
term in Euler’s equation is probably the time-derivative 
of the electric field. For ideal MHD, this term can be 
rewritten by first expressing E* in terms of the magnetic 
fields using the ideal MHD relation (4-8), and then us¬ 
ing (4-11) to eliminate the time-derivative of (see Zhang 
1989). This term is likely to be small in most applications; 
for example, it is 0{v‘^/c^) times smaller than the last two 
terms on the right hand side of (6-18). In such cases, ex¬ 
trapolating and iterating, or some other simple treatment, 
may be adequate to account for its contribution. 

It is instructive to take the Newtonian limit of equa¬ 
tion (6-18) and recover a familiar expression. With goo ^ 
-(1-1- 2(j)), where (p is the Newtonian potential, we find 

= pdiP. (6-19) 

In cartesian coordinates ( 7 ^^^ = 1), the Newtonian limit 

of equation (6-18) then becomes 

dtSi + dj{SiV^) = -diP - pdi4> + peEi 

-±d.{B^B,) + ^B^d,B, ( 6 - 20 ) 

or equivalently 

P^ = -V(F + PM)-/oV<^+ 7 ^-(B-V)B + peE, ( 6 - 21 ) 

at Att 

where we have defined the magnetic pressure 

Pm = ^ ( 6 - 22 ) 

87r 

and where V is the spatial gradient operator. Note that 
for a neutral plasma Pe = 0 the electric field P“ disappears 
entirely from the above Newtonian equation. 

7. SOURCE TERMS FOR THE GRAVITATIONAL FIELD 
EQUATIONS 

We now catalogue the source terms p (2-10), Si (2-11), 
Sij (2-12) and S (2-13) that appear in the Hamiltonian 
constraint (2-6), the momentum constraint (2-7) and the 
evolution equation (2-8). Inserting the fluid stress-energy 
tensor (5-1) into equations (2-10) - (2-13) yields the fluid 
contributions to the source terms: 


Pfluid — PqIiW P 

(7-1) 

= pohWu^ 

(7-2) 

cfluid cfluid 


cfluid _ p„,. , i 

+ pohW^ 

(7-3) 

Cfluid = 3P + po/i(VF2-1). 

(7-4) 


Next we assemble the electromagnetic contributions to 
the source terms. To do so, we first need to construct the 
electromagnetic stress-energy tensor from the Faraday 
tensor P“^, 

(7-5) 

Inserting (3-1), we first find 

FabF^’^ = 2(P,P* - P,P*), (7-6) 

where we have used eabc<^^^'^ = 2,'-^/. With = 

^bd^ce _ ^be^cd^ jjj ( 7 - 5 ) becomes 

F^^F\ = rfn'’E,E^ + 2n(“e^)=‘'PcSd 

- P“P'’ -h 7“*'P*PL (7-7) 
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Combining the last two equations then yields the electro¬ 
magnetic stress-energy tensor in 3 -I- 1 form 

47rTet = 

-b (7-8) 

This stress-energy tensor can now be inserted into (2- 
10) to (2-13) to obtain the electromagnetic source terms. 
For the mass-energy density pem we find 

dTrpem = = i(E,E* + B,B^) = ^(E^ + B^), 

(7-9) 

which is the energy density of the electromagnetic fields. 
The energy flux S'®™ reduces to the Poynting vector 

dTTSr = -TianbdTrT^t = 

= eijkE^B’^ = (E X B),. (7-10) 

The stress tensor Sf™ is 

dTTS®™ = = -E,E, - B,B, + + B2). 

(7-11) 

Its trace (2-13), finally, is equal to the mass-energy density 
Pem. 

4^Se„, = i(E2+B2). (7-12) 

The above results are not surprising: expressed in terms 
of the electromagnetic field components as measured by a 
normal observer, n“, i.e. an observer who is at rest with 
respect to the slices E, the 3-1-1 source terms have the 
same form as in flat space (compare exercise 5.1 in Mis- 
ner, Thorne & Wheeler 1973). 

8. COMPARISON WITH PREVIOUS TREATMENTS 

In this section we compare our notation and findings 
with those of Sloan & Smarr (1985, SS), Hawley & Evans 
(1988, EH; and 1989, HE) and Zhang (1989, Z). 

SS define the three-velocity ^SS by writing the four- 
velocity as 

=a-iW(l,a4s-/30 ( 8 - 1 ) 

and 

Ua = W{-a + vf) (8-2) 

(see equations (SS-2.1)). Here W is the Lorentz factor 
between and n“ 

W = —Uau'^ = au* (8-3) 

as in (4-3). The normalization u°'Ua = — 1 leads to 

W = (1 - (8-4) 

which shows that ^SS is the velocity of the fluid with re¬ 
spect to a normal observer. 

Zhang (1989) adopts the same formulation as SS, but 
denotes W with 7 (see equation (Z-2.11)). 

HE adopt the same definition of three-velocity as we 
do, defining the three-velocity u* to be the velocity with 
respect to coordinate observers, 

f(v = (8-5) 

(see equation (4-5) above). We use the subscript W since 
this definition is used in Wilson’s equations of relativistic 


hydrodynamics (see Wilson 1972; Hawley, Smarr & Wilson 
1984). With (8-5), the four-velocity can be written 

u“ = a-iW(l,i;V). ( 8 - 6 ) 

Comparing (8-1) and ( 8 - 6 ) shows that the two definitions 
of w* are related by 

~ (8-7) 

We can now compare the ideal MHD equation (4-8) in 

the different treatments. Since HE adopt the same defini¬ 

tion for u* = as we do, their equations (EH-A14) and 
(HE-14) should be identical to our equation (4-8). In their 
expression, however, the shift term is absent. This absent 
shift term can be traced back to their equation (HE-13), 
which does not agree with our equation (4-4). It is likely 
that the shift term was missed by dropping the term euj ■ 
The alignment of indices is incorrect in SS’s equation (2.9) 
(which they express in terms of m* instead of w*). Fixing it 
and utilizing (8-2) and (8-7) makes their equation equiva¬ 
lent to (4-8). 

To compare with Zhang’s ideal MHD equation (Z- 
2.12), we insert (8-7) into (4-8), which immediately yields 
Zhang’s result 

El = -eijkV^^B^, (8-8) 

showing that our result agrees with that of Zhang. 

We find similar errors in the Faraday equation. We 
found that the shift terms in (4-8) cancel all other shift 
terms when inserted into (3-10), ultimately yielding ex¬ 
pressions (4-11) and (4-13) which do not include any shift 
terms. With the shift terms being absent in equation (HE- 
14), the corresponding terms do not cancel, leading to the 
incorrect equations (HE-17) and (HE-18) (see also (EH- 
2.8) and (EH-A17)). 

Zhang’s expression for the Faraday equation (Z-2.13’) 
can be recovered by inserting (8-7) into (4-13) 

dtB^ = d, ((a4s - )B') , (8-9) 

which can be rewritten as 

= D, ((au^s - P^)B^ - {avis " 

^ ( 8 - 10 ) 

or 

= V X ((avss - /?) X B) . (8-11) 

This shows that our equations (4-11) and (4-13) again 
agree with the expressions of Zhang. 

Interestingly, Zhang refers to EH, and in fact their equa¬ 
tions look quite similar in that they both contain the above 
shift terms. However, Zhang uses ^SS as the three-velocity, 
while Hawley and Evans use Therefore the shift terms 
are correct in the former, but incorrect in the latter. 

Finally, we show that our equation (4-4) is equivalent to 
Zhang’s equation (Z-2.10). On the left hand side of (4-4) 
we rewrite 

Pe = UaJ^" = Ua{n°‘Pe + J“) = -WPe E U^E 

= H/(vss ■ J - Pe) (8-12) 

where we have used the decomposition (3-4). Inserting 
this and (8-7) into (4-4) yields 

J -b VE^(vss • j - Pe)vss = crlT(E -b vss X B), (8-13) 

which is identical to (Z-2.10). 
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9. SUMMARY 

We have assembled a complete set of Maxwell-Einstein- 
MHD equations, describing the structure and evolution of 
a relativistic, ideal MHD gas in a dynamical spacetime. 
We compare with previous treatments, and correct some 
errors in the existing literature. 

Our compilation of these equations is motivated by a 
large number of problems in relativistic astrophysics in 
which magnetic helds are likely to play an important role 
(see the incomplete list in Section 1). Self-consistent so¬ 
lutions to the Maxwell-Einstein-MHD equations will be 
necessary for a thorough understanding of these problems. 


and we therefore anticipate that relativistic MHD in dy¬ 
namical spacetimes will attract much interest in the future. 
We hope that our compilation of these equations will be 
useful for such investigations, particularly for treatments 
that will rely on numerical simulations. In Paper H, we 
use these equations to model the collapse of a magnetized 
star to a black hole. 
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5-I078I at the University of Illinois at Urbana-Champaign 
and NSF Grant PHY 0139907 at Bowdoin Gollege. 
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